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ABSTRACT 

The distribution function describing the acceleration of relativistic particles in an 
advection-dominated accretion disk is analyzed using a transport formalism that in- 
cludes first-order Fermi acceleration, advection, spatial diffusion, and the escape of par- 
ticles through the upper and lower surfaces of the disk. When a centrifugally-supported 
shock is present in the disk, the concentrated particle acceleration occurring in the 
vicinity of the shock channels a significant fraction of the binding energy of the accret- 
ing gas into a population of relativistic particles. These high-energy particles diffuse 
vertically through the disk and escape, carrying away both energy and entropy and al- 
lowing the remaining gas to accrete. The dynamical structure of the disk/shock system 
is computed self-consistently using a model previously developed by the authors that 
successfully accounts for the production of the observed relativistic outflows (jets) in 
M87 and Sgr A*. This ensures that the rate at which energy is carried away from the 
disk by the escaping relativistic particles is equal to the drop in the radial energy flux 
at the shock location, as required for energy conservation. We investigate the influ- 
ence of advection, diffusion, and acceleration on the particle distribution by computing 
the nonthermal Green's function, which displays a relatively flat power-law tail at high 
energies. We also obtain the energy distribution for the particles escaping from the 
disk, and we conclude by discussing the spectrum of the observable secondary radiation 
produced by the escaping particles. 
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1. INTRODUCTION 

Accretion flows onto supermassive black holes with mass M > 10 8 M & are believed to power 
the high-energy emission observed from quasars and active galactic nuclei (AGNs), as discussed 
by Lynden-Bell (1969; see also Novikov k Thorne 1973; Rees 1984; and Ford et al. 1994). When 
the accretion rate is relatively low compared with the Eddington value, the flow is radiatively 
inefficient, and the gas temperature approaches the virial value. In this case the disk has an 
advection-dominated structure, and most of the binding energy is swallowed by the black hole (e.g., 
Narayan, Kato, k Honma 1997; Blandford k Begelman 1999; Becker, Subramanian, k Kazanas 
2001; Becker k Le 2003). Despite the inefficiency of X-ray production in these sources, advection- 
dominated disks are luminous radio and 7-ray emitters, and they frequently display powerful bipolar 
jets of matter escaping from the central mass, presumably containing high-energy particles (e.g., 
Sambruna et al. 2004; Di Matteo et al. 2000; Allen, Di Matteo k Fabian 2000; Urry k Padovani 
1995; Owen, Eilek, k Kassim 2000). 

Although the effect of a standing shock in heating the gas in the post-shock region has been 
examined by a number of previous authors for both viscid (Chakrabarti k Das 2004; Lu, Gu, 
k Yuan 1999; Chakrabarti 1990) and inviscid (e.g., Lu k Yuan 1997, 1998; Yang k Kafatos 
1995; Abramowicz k Chakrabarti 1990) disks, the implications of the shock for the acceleration 
of nonthermal particles in the disk have not been considered in detail before. However, a great 
deal of attention has been focused on particle acceleration in the vicinity of supernova-driven 
shock waves as a possible explanation for the observed cosmic-ray energy spectrum (Blandford 
k Ostriker 1978; Jones k Ellison 1991). In the present paper we consider the analogous process 
occurring in hot, advection-dominated accretion flows (ADAFs) around black holes. These disks 
are ideal sites for first-order Fermi acceleration at shocks because the plasma is collisionless, and 
therefore a small fraction of the particles can gain a great deal of energy by repeatedly crossing the 
shock. Shock acceleration in the disk therefore provides an intriguing possible explanation for the 
powerful outflows of relativistic particles observed in many radio- loud systems (Le k Becker 2004) . 

The idea of shock acceleration in the environment of AGNs was first suggested by Blandford k 
Ostriker (1978). Subsequently, Protheroe k Kazanas (1983) and Kazanas k Ellison (1986) inves- 
tigated shocks in spherically-symmetric accretion flows as a possible explanation for the energetic 
radio and 7-radiation emitted by many AGNs. However, in these papers the acceleration of the 
particles was studied without the benefit of a detailed transport equation, and the assumption of 
spherical symmetry precluded the treatment of acceleration in disks. The state of the theory was 
advanced significantly by Webb k Bogdan (1987) and Spruit (1987), who employed a transport 
equation to solve for the distribution of energetic particles in a spherical accretion flow characterized 
by a self-similar velocity profile terminating at a standing shock. While more quantitative in nature 
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than the earlier models, these solutions are not directly applicable to disks since the geometry is 
spherical and the velocity distribution is inappropriate. Hence none of these previous models can 
be used to develop a single, global, self-consistent picture for the acceleration of relativistic particles 
in an accretion disk containing a shock. 

Accretion disks around black holes are not spherical, except perhaps in the innermost region, 
and therefore a transport equation written in cylindrical geometry is required in order to describe 
the acceleration of energetic particles in the accreting plasma. In this paper we develop a self- 
consistent cylindrically-symmetric model describing the acceleration of protons and/or electrons 
via first-order Fermi acceleration processes operating in a hot, advection-dominated accretion disk 
containing an isothermal shock. This is the third and final paper in a series studying the production 
of relativistic particles in ADAF disks with shocks. In Le & Becker (2004; hereafter Paper 1), 
we successfully applied our inflow/outflow model to explain the observed kinetic power of the 
outflows in M87 and Sgr A*, and in Le & Becker (2005; hereafter Paper 2), we discussed the 
detailed dynamics of ADAF disks containing isothermal shocks. Based on this dynamical model 
we can self-consistently compute the structure of the shocked disk and also the rate of escape of 
particles and energy from the disk at the shock location. In this paper we provide further details 
regarding the distribution in space and energy of the accelerated particles in the disk, and we also 
describe the energy distribution of the relativistic particles escaping from the disk at the shock 
location. Our focus here is on the computation of the Green's function for monoenergectic particle 
injection, which describes in detail how the particles are advected, diffused, and accelerated within 
the disk, resulting in a nonthermal particle distribution with a high-energy power-law tail. We 
also discuss the spectrum of the observable secondary radiation produced by the escaping particles, 
including radio (synchrotron) as well as inverse-Compton X-ray and 7-ray emission, which provides 
a quantitative basis for developing observational tests of the model. 

The first-order Fermi acceleration scenario that we focus on here (and in Papers 1 and 2) 
parallels the early studies of cosmic-ray acceleration in supernova shock waves. In analogy with the 
supernova case, we expect that shock acceleration in the disk will naturally produce a power-law 
energy distribution for the accelerated particles. Our earlier work in Paper 2 established that the 
first-order Fermi mechanism provides a very efficient means for accelerating the jet particles. In 
this paper, the solution for the Green's function, f G (E,r), describing the evolution of monoener- 
getic seed particles injected into the disk, is obtained by solving the transport equation using the 
method of separation of variables. The eigenvalues and the corresponding spatial eigenfunctions are 
determined using a numerical approach because the differential equation obtained for the spatial 
function cannot be solved in closed form. The resulting solution for the Green's function provides 
useful physical insight into how the relativistic particles are distributed in space and energy. We 
can also use the solution for f G (E,r) to check the self-consistency of the theory by comparing the 
numerical values obtained for the number and energy densities by integrating the Green's function 
with those obtained in Paper 2 by directly solving the associated differential equations for these 
quantities. 
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The remainder of the paper is organized as follows. In § 2 we briefly review the inviscid 
dynamical model used to describe the structure of the ADAF disk/shock system. In § 3 we develop 
and solve the associated steady-state particle transport equation to obtain the solution for the 
Green's function, f G (E,r), describing the evolution of monoenergetic seed particles injected from 
a source located at the shock radius, and in § 4 we evaluate the Green's function using dynamical 
parameters appropriate for M87 and Sgr A*. In § 5 we discuss the observational implications of 
our results for the production of secondary radiation by the outflowing relativistic particles, and in 
§ 6 we review our main conclusions. 

2. ACCRETION DYNAMICS 

It has been known for some time that inviscid accretion disks can display both shocked and 
shock- free (i.e., smooth) solutions depending on the values of the energy and angular momentum 
per unit mass in the gas supplied at a large radius (e.g., Chakrabarti 1989; Chakrabarti & Molteni 
1993; Kafatos & Yang 1994; Lu & Yuan 1997; Das, Chattopadhyay, & Chakrabarti 2001). Shocks 
can also exist in viscous disks if the viscosity is relatively low (Chakrabarti 1996; Lu, Gu, & Yuan 
1999), although smooth solutions are always possible for the same set of upstream parameters 
(Narayan, Kato, & Honma 1997; Chen, Abramowicz, & Lasota 1997). Hawley, Smarr, & Wilson 
(1984a, 1984b) have shown through general relativistic simulations that if the gas is falling with 
some rotation, then the centrifugal force can act as a "wall," triggering the formation of a shock. 
Furthermore, the possibility that shock instabilities may generate the quasi-periodic oscillations 
(QPOs) observed in some sources containing black holes has been pointed out by Chakrabarti, 
Acharyya, & Molteni (2004), Lanzafame, Molteni, & Chakrabarti (1998), Molteni, Sponholz, & 
Chakrabarti (1996), and Chakrabarti & Molteni (1995). Nevertheless, shocks are "optional" even 
when they are allowed, and one is always free to construct models that avoid them. However, in 
general the shock solution possesses a higher entropy content than the shock-free solution, and 
therefore we argue based on the second law of thermodynamics that when possible, the shocked 
solution represents the preferred mode of accretion (Becker & Kazanas 2001; Chakrabarti & Molteni 
1993). 

The possible presence of shocks in hot ADAF disks with low viscosity motivated us to explore in 
Paper 2 the relationship between the dynamical structure of the disk/shock system and the possible 
acceleration of relativistic particles that can escape from the disk to power the outflows commonly 
observed from radio-loud accretion flows around both stellar-mass and supermassive black holes. 
The disk/shock/outflow model is depicted schematically in Figure HJ In this scenario, the gas is 
accelerated gravitationally toward the central mass, and experiences a shock transition due to an 
obstruction near the event horizon. The obstruction is provided by the "centrifugal barrier," which 
is located between the inner and outer sonic points. Particles from the high-energy tail of the 
background Maxwellian distribution are accelerated at the shock discontinuity via the first-order 
Fermi mechanism, resulting in the formation of a nonthermal, relativistic particle distribution in 
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Fig. 1. — Schematic representation of the disk/shock/outflow model. The filled circles in the disk 
represent the accelerated particles, and the open circles represent the MHD scattering centers which 
advect with the background flow velocity. The seed particles are injected at the shock location. 

the disk. The spatial transport of the energetic particles within the disk is a stochastic process 
based on a three-dimensional random walk through the accreting background gas. Consequently, 
some of the accelerated particles diffuse to the disk surface and become unbound, escaping through 
the upper and lower edges of the cylindrical shock to form the outflow, while others diffuse outward 
radially through the disk or advect across the event horizon into the black hole. 

2.1. Transport Rates 

Becker & Le (2003) and Becker &: Subramanian (2005) demonstrated that three integrals of 
the flow are conserved in viscous ADAF disks when outflows are not occurring. The three integrals 
correspond to the mass transport rate 

M = 4irrHpu, (1) 

the angular momentum transport rate 

j = Mr 2 n - G , (2) 

and the energy transport rate 

E = -gn + M^-vl+ l -u 2 + ^ + ^ , (3) 

where p is the mass density, is the angular velocity, Q is the torque, H is the disk half-thickness, 
u is the positive radial inflow speed, = r Q is the azimuthal velocity, U is the internal energy 
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density, P = (7 — 1) U is the pressure, and 

_ , -GM 

$(r) = (4) 

r-r s 

denotes the pseudo-Newtonian gravitational potential for a black hole with mass M and Schwarzschild 
radius r s = 2GM/c 2 (Paczyhski & Wiita 1980). Each of the various densities and velocities repre- 
sents a vertical average over the disk structure. We also assume that the ratio of specific heats, 7, 
maintains a constant value throughout the flow. Note that all of the transport rates M, J, and E 
are defined to be positive for inflow. 

The quantities M and J remain constant throughout the entire flow, and therefore they rep- 
resent the rates at which mass and angular momentum, respectively, enter the black hole. The 
energy transport rate E also remains constant, except at the location of an isothermal shock if one 
is present in the disk. The torque Q is related to the gradient of the angular velocity O via the 
usual expression (e.g., Frank et al. 2002) 

g = -4 7 rr 3 Hp^, (5) 
dr 

where v is the kinematic viscosity. The disk half-thickness H is given by the standard hydrostatic 
prescription 

H(r) = (6) 
where c s represents the adiabatic sound speed, defined by 



Cs(r) EE ^-J , (7) 

and Qk denotes the Keplerian angular velocity of matter in a circular orbit at radius r in the 
pseudo-Newtonian potential (eq. [1]), defined by 

n> (r)= °" i« (8 ) 

r(r — r s ) z r dr 



2.2. Inviscid Flow Equations 

Following the approach taken in Paper 2, we focus here on the inviscid case since that allows the 
disk structure to be computed analytically, while retaining reasonable agreement with the dynamics 
of actual disks with low viscosities. We utilize the set of physical conservation equations employed 
by Chakrabarti (1989) and Abramowicz & Chakrabarti (1990), who investigated the structure of 
one-dimensional, inviscid, steady-state, axisymmetric accretion flows. In the inviscid limit, the 
torque Q vanishes, and the angular momentum per unit mass transported through the disk is given 

b y 

£(r) = — = r 2 £l(r) = constant . (9) 
M 
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We can combine equations ([3]) and ([7]) to show that the energy transported per unit mass in the 
inviscid case can be written as 

E I 2 I £ 2 c 2 s 

e = ^ = -u 2 -\ ~ H — + <£> (10) 

M 2 2 r 2 7 - 1 v ; 

where we have also utilized the relation I = rv^. The resulting disk/shock model depends on three 
free parameters, namely the energy transport per unit mass e, the specific heats ratio 7, and the 
specific angular momentum £. The value of e will jump at the location of an isothermal shock if 
one exists in the disk, but the value of I remains constant throughout the flow. This implies that 
the specific angular momentum of the particles escaping through the upper and lower surfaces of 
the cylindrical shock must be equal to the average value of the specific angular momentum for the 
particles remaining in the disk, and therefore the outflow exerts no torque on the disk (Becker, 
Subramanian, & Kazanas 2001). 

The values of the energy transport parameter e on the upstream and downstream sides of the 
isothermal shock at r = r* are denoted by e_ and e+, respectively. The condition e_ > e + must be 
satisfied as a consequence of the loss of energy through the upper and lower surfaces of the disk at 
the shock location. In order to satisfy global energy conservation, the rate at which energy escapes 
from the gas as it crosses the shock must be equal to the jet luminosity, i>j e tj and therefore we have 

Ljet = —AE = —MAe oc ergs s _1 , (11) 

where A denotes the difference between quantities on the downstream and upstream sides of the 
shock, so that Ae = e+ — e_ < 0. The drop in e at the shock has the effect of altering the transonic 
structure of the flow in the post-shock region, compared with the shock-free configuration (see 
Paper 2). In the inviscid case, the flow is purely adiabatic except at the location of an isothermal 
shock (if one is present), and therefore the pressure is related to the density profile via 

P = D p\ (12) 

where Do is a parameter related to the specific entropy that remains constant except at the location 
of a shock. The process of constructing a dynamical model for a specific disk/shock system begins 
with the selection of values for the fundamental parameters e_, £, and 7. For given values of these 
three quantities, a unique flow structure can be obtained using the procedure discussed in §§ 3 and 
4 of Paper 2, which yields a numerical solution for the inflow speed u(r). Following Narayan, Kato, 
& Honma (1997), we shall assume an approximate equipartition between the gas and magnetic 
pressures, and therefore we set 7 = 1.5. The range of acceptable values for e_ and £ is constrained 
by the observations of a specific object, as discussed in § 7 of Paper 2. 



3. STEADY-STATE PARTICLE TRANSPORT 

The gas in ADAF disks is hot and tenuous, and therefore essentially collisionless. In this 
situation, the relativistic particles diffuse mainly via collisions with magnetohydrodynamical (MHD) 
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scattering centers (e.g., Alfven waves), which control both the spatial transport and the acceleration 
of the particles. The magnetic scattering centers are advected along with the background flow, 
and therefore the shock width is expected to be comparable to the magnetic coherence length, 
Amag • Particles from the high-energy tail of the background Maxwellian distribution are accelerated 
at the shock discontinuity via the first-order Fermi mechanism, resulting in the formation of a 
nonthermal, relativistic particle distribution in the disk. The probability of multiple shock crossings 
decreases exponentially with the number of crossings, and the mean energy of the particles increases 
exponentially with the number of crossings. This combination of factors naturally gives rise to a 
power-law energy distribution, which is a general characteristic of Fermi processes (Fermi 1954). 

Two effects limit the maximum energy that can be achieved by the particles. First, at very 
high energies the particles will tend to lose energy to the waves due to recoil. Second, the mean 
free path A mag will eventually exceed the thickness of the disk as the particle energy is increased, 
resulting in escape from the disk without further acceleration. As a consequence of the random 
walk, some of the accelerated particles diffuse to the disk surface and become unbound, escaping 
through the upper and lower edges of the cylindrical shock to form the outflow, while others diffuse 
outward radially through the disk or advect across the event horizon into the black hole (see Fig. [I]). 
Our goal in this paper is to analyze the Green's function f G (E, r) in the disk based on the transport 
equation and the dynamical results discussed in Paper 2. We focus specifically on models 2 and 5 
from Paper 2, which describe the disk structures for M87 and Sgr A*, respectively. 



3.1. Transport Equation 

The particle transport formalism we utilize includes advection, spatial diffusion, first-order 
Fermi acceleration, and particle escape. Our treatment of the Fermi process includes both the 
general compression related to the overall convergence of the accretion flow, as well as the enhanced 
compression that occurs at the shock. To avoid unnecessary mathematical complexity, we utilize 
a simplified model for the spatial transport in which only the radial (r) component is treated in 
detail. In this approach, the diffusion and escape of the particles in the vertical (z) direction is 
modeled using an escape-probability formalism. We will adopt the test-particle approximation 
utilized in Paper 2, meaning that the dynamical effect of the relativistic particle pressure on the 
flow structure is assumed to be negligible. This assumption is valid provided the pressure of the 
relativistic particles is much smaller than the background (thermal) pressure, which is verified for 
models 2 and 5 in § 8 of Paper 2. We also assume that the injection of the seed particles and the 
escape of the accelerated particles from the disk are both localized at the isothermal shock radius, 
which is consistent with the strong concentration of first-order Fermi acceleration in the vicinity of 
the shock. 

The localization of the particle injection and escape at the shock radius allows us to make a 
direct connection between the jump in the energy flux in the disk at the shock location and the 
energy carried away by the escaping particles. This connection is essential in order to maintain self- 
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consistency between the dynamical and transport calculations, as discussed in Paper 2. In a steady- 
state situation, the Green's function, f G (E,r), representing the particle distribution resulting from 
the continual injection of monoenergetic seed particles from a source located at the shock radius 
(r = r*) satisfies the transport equation (Becker 1992) 

9f G n ^ * Id 







V • F - ^ m ^V-Vf Q ) + f source - /esc , (13) 



dt 

where the specific flux F is evaluated using 

and the quantities E, k, and v represent the particle energy, the spatial diffusion coefficient, and 
the vector velocity, respectively. The velocity has components given by v = v r r + v z z + (ft, with 
v r = —u < since we are dealing with inflow. 

The injection of the monoenergetic seed particles is described by the function 

■ _ N S(E-E )5(r-n) 

/source- (4^)2^ > W 

where Eq is the energy of the injected particles, Nq is the particle injection rate, and H* = H(r*) 
represents the half-thickness of the disk at the shock location. The escape of the particles through 
the upper and lower surfaces of the disk is represented by the term 

/ esc = A c5(r - r*) / G , (16) 

where the dimensionless parameter is computed using (see Appendix B of Paper 2) 

3/£* 



and = (k- + k+)/2 denotes the mean value of the spatial diffusion coefficient at the shock 
radius. The subscripts "-" and "+" will be used to denote quantities measured just upstream and 
just downstream from the shock, respectively. The condition Aq < 1 is required for the validity of 
the diffusive treatment of the vertical escape employed in our approach. 

In order to compute the Green's function f G (E,r), we must specify the injection energy of 
the seed particles Eq as well as their injection rate Nq. Following the approach taken in Paper 2, 
we set the injection energy using Eq = 0.002 ergs, which corresponds to an injected Lorentz factor 
Tq = EQ/{mp(?) ~ 1.3, where m p is the proton mass. Particles injected with energy Eq are 
subsequently accelerated to much higher energies due to repeated shock crossings. We find that 
the speed of the injected particles, vq = c(l — Tq 2 ) 1 / 2 , is several times larger than the mean ion 
thermal velocity at the shock location, v rms = (3/cT^/mp) 1 / 2 , where T* is the ion temperature at the 
shock. The seed particles are therefore picked up from the high-energy tail of the local Maxwellian 
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distribution. With Eq specified, we can compute the particle injection rate iVo using the energy 
conservation condition 

N E = -MAe ocergss -1 , (18) 

where Ae < represents the jump in the energy transport per unit mass as the plasma crosses 
the shock (see eqs. [10] and [H]). This relation ensures that the energy injection rate for the seed 
particles is equal to the energy loss rate for the background gas at the isothermal shock location. 

The total number and energy densities of the relativistic particles, denoted by n r and U r , 
respectively, are related to the Green's function via 

/■oo roc 

n r (r) = AnE 2 f G (E, r) dE , U r (r) = 4ttE 3 f G (E,r) dE , (19) 
io io 

which determine the normalization of f G . Equations (|13|) and (|14p can be combined to obtain the 
alternative form 

^•V/ G = ||^V-zT+V - («V/ G ) + / source - Asc , (20) 

where the left-hand side represents the co-moving (advective) time derivative and the terms on the 
right-hand side describe first-order Fermi acceleration, spatial diffusion, the particle source, and the 
escape of particles from the disk at the shock location, respectively. Note that escape and particle 
injection are localized to the shock radius due to the presence of the 5-functions in equations (fT5j) 
and (|16p . Our focus here is on the first-order Fermi acceleration of relativistic particles at a 
standing shock in an accretion disk, and therefore equation ([20]) does not include second-order 
Fermi processes that may also occur in the flow due to MHD turbulence (e.g., Schlickeiser 1989a,b; 
Subramanian, Becker, &: Kazanas 1999; Becker et al. 2006). Since no energy loss mechanisms are 
included in the model, all of the particles injected with energy Eq will experience an increase in 
energy. 



Under the assumption of cylindrical symmetry, equations () 1 5 j) . (|16p . and (|20p can be combined 
to obtain 



df G df G E df G 



l_d_ , , | dv z 

r dr r dz 



dr dz 3 dE 

N 8{E - Eq) 5{r - r* 



1 d 



rn 



{4nE ) 2 nH* 



r dr \ dr 
AQ C 8{r-r*)f G , (21) 



where the escape of particles from the disk is described by the final term on the right-hand side. 
We focus here on the vertically-integrated form of the transport equation, which can be written as 
(see Paper 2, Appendix A) 

df G 1 d , „ . Edf G t 1 / df G 



dr r dr^ ^3 dE r dr \ dr 
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where u = —v r and / G and k are vertically averaged quantities. Due to the presence of the velocity 
derivative on the right-hand side of equation (|22p , the first-order Fermi acceleration of the particles 
is most pronounced near the shock, where u is discontinuous. In the vicinity of the shock, we find 
that 

du 

> (u- - u + ) 8(r - r*) , r -> r* , (23) 

ar 

where U- and u+ denote the positive inflow speeds just upstream and downstream from the shock, 
respectively. The velocity jump condition appropriate for an isothermal shock is (see Paper 2) 

£ = ^ w» 

where A4- is the incident Mach number on the upstream side of the shock. Although the Fermi 
acceleration of the particles is concentrated at the shock, the rest of the disk also contributes to the 
particle acceleration because of the general convergence of the MHD waves in the accretion flow. 

In order to close the system of equations and solve for the Green's function using equation (I22p . 
we must also specify the radial variation of the diffusion coefficient k. The behavior of k can be 
constrained by considering the fundamental physical principles governing accretion onto a black 
hole, which we have demonstrated in Paper 2. The precise functional form for the spatial variation 
of k is not completely understood in the accretion disk environment. In order to obtain a mathe- 
matically tractable set of equations with a reasonable physical behavior, we utilize the general form 
(see Paper 2, § 5) 

«(r) = k u(r) r s ( — - 1 ) , (25) 



where kq > is a dimensionless constant that can be computed for a given source based on energy 
considerations (see § 4) . Due to the appearance of the radial speed u in equation (j25j) , we note that 
K exhibits a jump at the shock. This is expected on physical grounds since the MHD waves that 
scatter the ions are swept along with the thermal background flow, and therefore they should also 
experience a density compression at the shock. We also point out that the dimensionless escape 
parameter Aq can be evaluated by combining equations (|17j) and (|25|) . which yields 

*=(^)'(H 4<1 ' (26> 

where u* = (u- + n+)/2 denotes the mean value of the inflow speed at the shock location. Note 
that the value of the diffusion parameter kq is constrained by the inequality in equation (j26|) . 



3.2. Separation of Variables 

For values of the particle energy E > Eq, the source term in equation (I22j) vanishes and the 
remaining equation is separable in energy and space using the functions 

f n (E,r)=(^\ Xn Y n (r), (27) 
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where A n are the eigenvalues, and the spatial eigenfunctions Y n (r) satisfy the second-order ordinary 
differential equation 



Hu d^ = ^d y Id 

dr 6r dr r dr 



rHn ■ 



dr 



A cH*5(r - r*)Y n 



(28) 



The eigenvalues \ n are determined by applying suitable boundary conditions to the spatial eigen- 
functions, as discussed below. Once a numerical solution for the inflow speed u(r) has been obtained 
following the procedures described in Paper 2, we can compute H(r) and k(t) using equations ([6]) 
and (|25p . respectively. Since the coefficients in equation (128j) cannot be represented in closed form, 
analytical solutions cannot be obtained and therefore the eigenfunctions Y n (r) must be determined 
numerically. Away from the shock (r / r t ), equation (|28|) reduces to 



d 2 Y n 
dr 2 



+ 



d\a{rHK) u 
dr k 



dY n X n uY n dln(rHu) 



dr 



+ 



3 K 



which can be rewritten as 



d 2 Y n 
dr 2 



+ 



dln(rHu) 
dr 



+ 



+ 



K (r 



dYn 
dr 



+ 



dr 



^n^g Yn 

3n (r - r s y 







dln(rHu) 
dr 



, 



(29) 



(30) 



where we have used equation (|25p to substitute for k. Once a solution for the inflow speed u(r) 
has been determined using the procedure described in § 2, the variation of the disk half-thickness 
H(r) can be computed by combining equations (H|), ©, and (fT2|) . The result obtained is 



H(r) 



1 



Q K (r) 



,3/2/ 



i u(r) 



K 



(l-7)/(l+7) 



(31) 



where is given by equation ([8]) and K is the "entropy parameter", which remains constant 
except at the location of an isothermal shock if one is present in the flow (see Paper 2 and Becker 
&; Le 2003). By utilizing the solutions obtained for u{r) and H(r), we are able to compute all of 
the coefficients in equation (|30|) . which allows us to solve numerically for the spatial eigenfunctions 
Y n (r). 



3.3. Eigenvalues and Eigenfunctions 



The global solution for the eigenfunction Y n (r) must satisfy continuity and derivative jump 
conditions associated with the presence of the shock/source at radius r = r*. By integrating 
equation (|28p with respect to radius in the vicinity of the shock, we find that 



A 



— uY n + K 
6 dr 



A[y n 

dY r 







-cA)^n(r*) , 



(32) 
(33) 
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where the symbol A denotes the difference between postshock and preshock quantities, and we 
have made use of the fact that H is continuous across the isothermal shock. Equations (|32p and 
(|33p establish that Y n (r) must be continuous at the shock location, and its derivative must display 
a jump there. 

Since the spatial eigenfunctions Y n (r) satisfy the second-order linear differential equation (I30p . 
we must also impose two boundary conditions in order to determine the global solutions and the 
associated eigenvalues A n . The required conditions can be developed based on our knowledge of 
the dynamical structure of the disk near the event horizon (r = r s ) and at large radii (see Becker 
& Le 2003 and § 4 of Paper 2). By utilizing the Frobenius expansion method to identify the 
dominant terms in equation (|30p . we find that the global solutions for the spatial eigenfunctions 
can be written as 

' G™(r) , r < r* , 



Y n (r) 



(34) 

a n G° n nt (r) , r > r* , 



where a n is a matching constant and G™(r) and G° ut (r) denote the fundamental solutions to 
equation (j30j) possessing the asymptotic behaviors 

, -A n /(3 7 +3) 



(35) 



at small and large radii, respectively. The value of a n is determined using the continuity condition 
(see eq. [32]) , which yields 

an = G^W) ■ (36) 

The validity of the asymptotic forms in equations (|35p is verified in § 4.1 via comparison with the 
numerical solutions obtained for the spatial eigenfunctions. 

We use a bi-directional integration scheme to solve for the eigenvalues A n based on the bound- 
ary conditions in equations ([35]) . The value of A n is iterated using a bisection method until the 
Wronskian of the inner and outer solutions vanishes at a matching radius located in the postshock 
region. Once a particular eigenvalue has been determined to the required accuracy, the matching 
constant a n is computed using equation (f36l) . This procedure is repeated until the desired number 
of eigenvalues and eigenfunctions have been obtained. The sequences of eigenvalues associated with 
the model 2 (M87) and model 5 (Sgr A*) parameters are plotted in Figure [2j Note that Ai ~ 4 in 
all cases, indicating that the acceleration is quite efficient, in analogy with the case of cosmic-ray 
acceleration (e.g., Blandford & Ostriker 1978). The first two eigenvalues in models 2 and 5 are 
Ai = 4.165, A2 = 6.415 and Ai = 4.180, A2 = 6.344, respectively. It follows that at high energies 
the power-law slope of the energy distribution is dominated by the first eigenvalue Ai, since all of 
the subsequent eigenvalues are much larger (see Fig. [2]). 
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Fig. 2. — Eigenvalues for model 2 (M87; filled squares) and model 5 (Sgr A*; filled circles). 



3.4. Eigenfunction Orthogonality 

We can establish several useful general properties of the eigenfunctions by rewriting equa- 
tion (1301) in the Sturm-Liouville form 



d 
dr 



S(r, 



where 



dYn 

dr 
{ 1 



+ \ n uj(r)Y n (r) = , 



_ rHn 
S{r) = — exp< — 



-i 



-i 



and the weight function tj(r) is defined by 

u;(r) e 



uS" d\n{rHu) 



(37) 
(38) 

(39) 



3 k dr 

Note that tcj(r) displays a <5- function discontinuity at the shock through its dependence on the 
derivative of u(r). In the vicinity of the shock, we can combine equations (I23p . (|38p . and (j39H to 
show that 



io(r) 



5(r - r*) , 



(40) 



3k^ J 

Based on the boundary conditions satisfied by the spatial eigenfunctions Y n {r) (see eqs. [35]), we 

demonstrate in Appendix A that the eigenfunctions satisfy the orthogonality relation 

poo 

Y n (r)Y m (r)uj(r)dr = , m ^ n . (41) 

' s 
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3.5. Eigenfunction Expansion 

We have established that Y n (r) is the solution to a standard Sturm-Liouville problem, and 
therefore it follows that the set of eigenfunctions is complete. Consequently, the Green's function 
f G (E,r) can be represented using the eigenfunction expansion (see eq. [27] ) 

f G (E,r)=J2b n Y n (r)i^—j , E > E , (42) 

where b n are the expansion coefficients (with either positive or negative signs) and the value of 
Af max is established by analyzing the term-by-term convergence of the series. We utilize the value 
Af max = 10 in our numerical examples, which generally yields an accuracy of at least three decimal 
digits. Note that f G (E,r) = for all E < Eq because no deceleration processes are included in the 
particle transport model. 

The expansion coefficients bi, 62,63, . . . can be calculated by utilizing the orthogonality of the 
spatial eigenfunctions. At the source energy, E = Eq, equation (|42p reduces to 

-/Vmax 

f G (E ,r)= 2 b m Y m (r) . (43) 

m=l 

Multiplying each side of this expression by the product Y n {r) u){r) and integrating from r = r s to 
r = 00 yields 

/OO ^Vmax /*00 

f G (E ,r)Y n (r)u(r)dr = ^ b m / Y m (r) Y n {r)u(r)dr . (44) 
s m=l Jr s 

Based on the orthogonality of the eigenfunctions, as expressed by equation (|41|) . we note that only 
the m = n term survives on the right-hand side of equation (|44p . and we therefore obtain 

poo 

f G (E ,r)Y n (r)uj(r)dr = b n Y^{r)uj(r)dr . (45) 



G 1 

S 



Solving this expression for the expansion coefficient b n yields 

J r °°f G (E ,r)Y n (r)u;(r)dr 

bn = — = 

where the quadratic normalization integrals, l n , are defined by 



(46) 



nOO 

J„ = / Y*(r)u(r)dr . (47) 

s 



To complete the calculation of the expansion coefficients, we need to evaluate the distribution 
function at the source energy, f G (Eo,r). This can be accomplished by using equation ([23]) to 
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substitute for the velocity derivative in equation (|22p and then integrating with respect to E in a 
small range around the injection energy Eq, which yields 



3iV 



{4m) 2 Elr*H*(u--u + ) ' ' (48) 

. , r ± r* . 



/ G (^o,r) 



Although the source term in the transport equation (|2ip is a delta function, we note that the 
Green's function evaluated at the injection energy Eq is not a delta function. Instead, as indicated 
by equation (148 h . we find that the Green's function at the injection energy is zero for radii away 
from the shock, and it is equal to a finite value at the shock (r = r*). From a physical point of 
view, this reflects the fact that the particle injection occurs at the shock location in our model. 
The strong acceleration experienced by the particles at the shock eliminates the appearance of a 
delta-function behavior in the resulting Green's function. It is interesting to compare this with 
the behavior observed in the case of particle acceleration in a plane-parallel supernova shock wave, 
which was analyzed in detail by Blandford & Ostriker (1978). Our results are similar to theirs, 
except that the Blandford & Ostriker Green's function evaluated at the injection energy Eq has a 
finite value at all locations in the flow, because there is no convergence in the plasma surrounding 
the plane-parallel supernova shock wave. However, in our model, the plasma converges at all radii 
in the accretion flow (not just at the shock), and this additional particle acceleration causes the 
Green's function to vanish at the injection energy for all r ^ r*. 

Utilizing equation ([4*8]) to substitute for f G (Eo,r) in equation (j4"6]) and carrying out the inte- 
gration, we obtain 

K = , , (49) 

where we have made use of the fact that the weight function u{r) has a 5- function behavior close to 
the shock (see eq. [30] ) . The singular nature of the weight function also needs to be considered when 
computing the normalization integrals I n defined in equation (|4T|) . By combining equations (f4TJj) 
and (j4"T|) . we find that 



r* — e 



l n = lim / u/(r) F» dr + / u(r) Y*(r) dr + ~ + ) Y^(r*) , (50) 
and we shall use this expression to evaluate the normalization integrals. 



4. ASTROPHYSICAL APPLICATIONS 

In Paper 2 we investigated the acceleration of particles in an inviscid ADAF disk containing 
an isothermal shock. For a given source with a known black hole mass M, accretion rate M, and 
jet kinetic power Lj et , we found that a family of flow solutions can be obtained, with each solution 
corresponding to a different value of the diffusion parameter kq (see eq. [25] ) . We concluded that 
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Fig. 3.— The fundamental solutions to eq. ([28]) G^r) and G^ ut (r) (so/id Zmes) for model 2 (M87) 
are compared with the corresponding asymptotic solutions <?i n (r) and g° ut (V) (dashed lines) in panels 
a and 6, respectively (see eqs. [32] )• The associated global solution for the first eigenfunction Y\(r) 
(eq. [3l]) is plotted in panel c. The shock location at r = r» is indicated. 



the upstream energy flux, e_, is maximized for a certain value of kq, and we adopted that particular 
value for kq in our models for M87 and Sgr A* in Paper 2. In Table Q] we list the values of the 
parameters M, M, Lj et , e_, £, e+, r*, and associated with models 2 and 5, which describe M87 
and Sgr A*, respectively. 



4.1. Numerical Solutions for the Eigenfunctions 

In order to compute the Green's function f G (E, r) , we must first solve for the eigenvalues A n and 
eigenfunctions Y n (r) following the procedures described in §§ 3.3 and 3.4 based on equation (|30p . 
Once the eigenvalues and eigenfunctions have been determined, we can compute the expansion 
coefficients b n using equation (09]), after which f G (E, r) is evaluated using equation ([4^2]) . In Figure[3] 
we compare the fundamental asymptotic solutions G™(r) and G° ut (r) with the asymptotic solutions 



Table 1. Disk Structure Parameters 



Source 


Model 


M 


M 

(M yr-1) 


(ergs s~ V) 


e_ 


£ 






C 


M87 
Sgr A* 


2 
5 


3.0 x 10 9 
2.6 x 10 6 


1.3 x lO" 1 
8.8 x 10~ 7 


5.5 x 10 43 
5.0 x 10 38 


1.527 x 10~ 3 
1.229 x 10~ 3 


3.1340 
3.1524 


-5.746 x 10" 3 
-8.749 x 10~ 3 


21.654 
15.583 


0.108 
0.124 



Note. - All quantities are expressed in gravitational units (GM = c = 1) except as indicated. 
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Fig. 4. — Same as Figure [3] but for model 5, which describes Sgr A* 
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gjf(r) and <7^ ut (r) for model 2 with n = 1 (see eqs. [35]). The corresponding results for model 5 
are plotted in Figure [H The excellent agreement between the G{r) and g{r) functions confirms 
the validity of the asymptotic relations we have employed near the event horizon and also at large 
radii. Figures [3] and H] also include plots of the global solutions for the first eigenfunction Y\{r) 
associated with models 2 and 5, respectively (see eq. [S]). Note that the global solutions display 
a derivative jump at the shock location, as required according to equation (j33|) . 



4.2. Green's Function Particle Distribution 

We can combine our results for the eigenvalues, eigenfunctions, and expansion coefficients to 
evaluate the Green's function f G (E,r) using equation (j4"2j) . In Figures [5^ and [5b we plot f G (E,r) 
as a function of the particle energy E at various radii r in the disk for models 2 and 5, respectively. 
Note that at the injection energy (E = Eq), the Green's function is equal to zero except at the 
source/shock radius, r = r*, in agreement with equation (08]). This reflects the fact that the 
particles are rapidly accelerated after they are injected, and therefore all particles have energy 
E > Eq once they propagate away from the shock location. The positive slopes observed at low 
energies when r ^ r* result from the summation of terms with either positive or negative expansion 
coefficients b n (see § 3.5). At energies above the turnover, the spectrum has a power-law shape 
determined by the first eigenvalue, Ai ~ 4 (see Fig. 2). The relatively flat slope of the Green's 
function above the turnover reflects the high efficiency of the particle acceleration in the shocked 
disk. By contrast, in situations involving weak acceleration, the Green's function has a strong peak 
at the injection energy surrounded by steep wings (e.g., Titarchuk & Zannias 1998). 

In Paper 2 we established that most of the injected particles are advected inward by the 
accreting MHD waves, eventually crossing the event horizon into the black hole. Hence only a 
small fraction of the particles are able to diffuse upstream to larger radii, which is confirmed by the 
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Fig. 5. — Results for the relativistic particle Green's function f G (E, r) in ergs -3 cm -3 computed 
using eq. (|42h for (a) M87 (model 2), and (b) Sgr A* (model 5). The value of the radius r in units 
of GM/c 2 is indicated for each curve, with r = denoting the shock location. 



strong attenuation of the particle spectrum with increasing r apparent in Figure El The particles in 
the inner region (r < r*) exhibit the greatest overall energy gain because they are able to cross the 
shock multiple times, and they also experience the strong compression of the flow near the event 
horizon. This is consistent with the results for the mean energy distribution, (E) = U r /n r , plotted 
in Figure 9 of Paper 2. 



4.3. Number and Energy Density Distributions 

In our numerical examples, we have generated results for the Green's function f G (E,r) by 
summing the first ten eigenfunctions using equation (142 j) . Once the Green's function energy dis- 
tribution has been determined at a given radius, we can integrate it with respect to the particle 
energy E to obtain the corresponding values for the number and energy densities. By combining 
equations (fl"9j) and ([4l?|) and integrating term-by-term, we obtain 

n?(r) s 4vr£ 3 £ , ^ (r) s 4^ £ . (51) 

n=l An 6 n=l An 4 

Equations (I5ip provide an important tool for checking the self-consistency of our entire formalism. 
This is accomplished by comparing the results for the number and energy densities computed using 
equations (|5ip with the corresponding values obtained by directly solving the differential equations 



-20- 




1 2 5 10 20 50 100 1 2 5 10 20 50 100 

r r 

Fig. 6. — Plots of (a) the relativistic particle number density and (b) the relativistic particle energy 
density in cgs units for model 2, which describes M87. The solid lines represent the solutions 
computed by numerically integrating the associated differential equations, and the filled circles 
represent the corresponding results obtained by integrating the Green's function using eqs. f)5 1 [) . 
The location of the shock at r = r t is indicated. Note the close agreement between the results, 
which confirms the accuracy of our computational model. 

for n r and U r using the method developed in Paper 2. In principle, the two sets of results should 
agree closely if our procedure for computing f G (E,r) is robust. We compare the model 2 results 
obtained for n r and U r in Paper 2 with those computed using equations (|5ip in Figure El The 
corresponding comparison for model 5 is carried out in Figure Note the close agreement between 
the two sets of results, which confirms the validity of the analysis involved in searching for the 
eigenvalues, computing the eigenfunctions, solving for the expansion coefficients, and calculating 
the Green's function. Moreover, the results in Figures 6 and 7 indicate that equation (|42p converges 
successfully using the first ten terms (iV max = 10) in the series. 



4.4. Escaping Particle Distribution 

The Green's function f G (E, r) computed using equation (|4"2"1) represents the energy distribution 
of the relativistic particles inside the disk at radius r. However, our primary goal in this paper is 
the calculation of the energy distribution of the relativistic particles escaping from the disk at the 
shock radius, r*. In order to compute the energy spectrum of the escaping particles, we need to 
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Fig. 7. — Same as Figure [6] but for model 5, which describes Sgr A*. 

integrate equation (I16p over the volume of the disk, which yields 

A>f c = (4TrE) 2 nH*cA f G (E,n) , (52) 

where AJ? C dE represents the number of particles escaping from the disk per unit time with energy 
between E and E + dE. At the highest particle energies, the escaping particle number spectrum 
has a power-law shape, with A|? c oc E~ a , where a = Ai — 2 and Ai is the first eigenvalue (see 
eqs. [12] and [52]). The total number of particles escaping from the disk per second, N esc , is related 
to the spectrum AJ? C via 

/•oo 

N esc = N^ c dE = 4irr*H*cA n* , (53) 
J o 

where n* = n r (r*) is the relativistic particle number density at the shock location. The correspond- 
ing result for the total energy escape rate is given by 

/•oo 

L csc = / iVf c EdE = AirnH^ cA U* , (54) 
Jo 

where U* = U r (r#) is the relativistic particle energy density at the shock location. 

The escaping particle distributions for M87 (model 2) and Sgr A* (model 5) are plotted in 
Figure [HJ The results obtained in the case of M87 are about five orders of magnitude larger than 
those associated with Sgr A* due to the corresponding difference in the reported jet luminosity 
Ljet for these two sources (see Table [1]). In Table [2] we list the values obtained for Ao, Ae SC , n*, 
£/*, k*, kq, Aq, E esc , and in models 2 and 5, where E'esc = U*/n* denotes the mean energy 
of the particles escaping from the disk to power the jet, and Too represents the terminal Lorenz 
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Fig. 8. — Plots of (a) the number distribution N^ c and (b) the energy distribution E iVJ? c for the 
relativistic particles escaping at the shock location, r = r* (see eq. [52]). The solid curves represent 
model 2 (M87), and the dash curves represent model 5 (Sgr A*). 



factor of the jet. We note that the values obtained for L esc using the integral in equation (|54p 
are in excellent agreement with the values for Lj e t reported in Paper 2 and listed in Table [TJ The 
agreement between L esc and Lj et confirms that our model satisfies global energy conservation. We 
also find that the values for iV esc obtained using the integral in equation (|53[) are in excellent 
agreement with the values listed in Table [2j which were taken from Paper 2. These tests confirm 
the validity of the analytical approach we have utilized in our derivation of the Green's function, 
given by equation (|42|) . 



Table 2. Transport Equation Parameters 



Model 


N 
(s-1) 


JV oac 
No 


n, 

(cm - ^) 


(ergs cm - ^) 


re* 




A 


-Bosc 

Eq 


Too 


2 
5 


2.75 x 10 46 
2.51 x 10 41 


0.17 
0.18 


2.01 x 10 4 
4.33 x 10 5 


2.39 x 10 2 
4.71 x 10 3 


0.427877 
0.321414 


0.02044 
0.02819 


0.0124 
0.0158 


5.95 
5.45 


7.92 
7.26 



Note. - All quantities are expressed in gravitational units (GM = c = 1) except as indicated. 
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5. OBSERVATIONAL IMPLICATIONS 

The jets of relativistic particles produced by the first-order Fermi shock acceleration mechanism 
considered here will flow outwards from the disk with initially high pressure, and low velocity. The 
rapid expansion of the plasma will lead to strong acceleration as the gas cools and the random 
internal energy is converted into ordered kinetic energy (see Paper 2). The terminal Lorenz factor, 
Too, associated with the M87 (model 2) and Sgr A* (model 5) jets considered here are found to 
be 7.92 and 7.26, respectively, based on comparison with the observed jet kinetic luminosities. 
However, robust confirmation of the model requires a detailed comparison between the predicted 
emission and the observed multifrequency spectra for a variety of sources. The jets envisioned 
here are expected to be charge neutral, with the protons "dragging" along an equal number of 
electrons. The electrons are expected to cool rapidly via inverse-Compton and bremsstrahlung 
emission produced while the particles are still close to the active nucleus. On the other hand, the 
relativistic protons will cool inefficiently via direct radiation, and they also lose very little energy 
to the electrons via Coulomb coupling due to the low gas density. Hence it is expected that the 
ions will maintain their high energies until they collide with particles or radiation in the ambient 
medium. Detailed consideration of this issue is beyond the scope of the present paper, and will 
be pursued in future work. However, we briefly discuss a few of the probable scenarios for the 
production of the observed radio, X-ray, and 7-ray emission by the jet particles below. 

Scenario 1: Anyakoha et al. (1987) and Morrison et al. (1984) have suggested that high-energy 
protons in the jets may undergo proton-proton collisions with the ambient medium, resulting in the 
production of charged and neutral pions. The neutral pions subsequently decay into 7-rays, and the 
charged pions decay into electron-positron pairs which produce further 7-rays via inverse-Compton 
and annihilation radiation. These processes can explain the observed radio and 7-ray emission in 
blazar sources such as 3C 273. Anyakoha et al. (1987) demonstrated that the production of the 
observed 7-rays requires protons with Lorentz factors in the range 1.4 < 7 < 86, and the production 
of the observed radio emission requires protons with Lorentz factors in the range 7 < 7 < 300. Based 
on our model, we find that the terminal proton Lorentz factor, F^, for the jets in M87 and Sgr A* is 
given by F^ ~ 8 and ~ 7, respectively (see Table [2|). We therefore conclude that the escaping 
protons accelerated in our model reach energies high enough to produce the observed radio and 
7-ray emission in M87 and Sgr A*. In addition to the radiation emitted by the escaping, outflowing 
protons, the relativistic protons remaining inside the disk will produce additional 7-rays (e.g., Eilek 
&: Kafatos 1983; Bhattacharyya et al. 2003, 2006), although the escape of these photons from the 
disk may be problematic due to 7-7 attenuation (Becker &; Kafatos 1995). 

Scenario 2: Usov & Smolsky (1998) suggested that the nonthermal 7-ray and X-ray emission 
from AGN jets may be produced by energetic electrons accelerated due to collisions between rela- 
tivistically moving jet clouds and the ambient medium. The collisions accelerate the electrons up to 
energies ~ m p c 2 F^, where is the terminal Lorentz factor of the jet in our model. Gamma-rays 
are produced when the high-energy electrons scatter external radiation. Based on the results pre- 
sented by Usov <fc Smolsky (1998), we find that the mean energy of the scattered photons emerging 
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from the jet in this scenario is <e 7 ) < 400r^,keV, where Too ~ 7 — 8 (see Table [2|). Additional 
X-rays and 7-rays may also be produced as a result of direct upscattering of background photons 
by beamed electrons if the Lorentz factor ~ 10, which is comparable to the range predicted by 
our model (e.g., Dermer & Schlickeiser 1993). 

Scenario 3: This last scenario is quite different from the first two. Sikora & Madejski (2000) 
argue that the jet in optically violently variable (OVV) quasars contains a mixture of compositions, 
since the pure e + e~ jets overproduce soft X-rays and the pure e — p jets underproduce nonthermal 
X-rays. To resolve these issues, they suggest that initially e — p plasma is injected at the base of the 
jets, with the protons providing the inertia to account for the kinetic luminosity of the jet. This 
is consistent with our model of proton/electron escape at the shock. Subsequently, e + e~ pairs 
are produced via interactions with the hard X-rays/soft 7-rays emitted from the hot accretion disk 
corona. This leads to velocity and density perturbations in the jet and also to the formation of 
shocks, where the particles are further accelerated. The consequence of the Sikora & Madejski 
(2000) model is that pair production in the accretion-disk corona can be responsible for the fast 
(~ day) variability observed in OVV quasars at very small angles relative to the jet axis, and also 
detected in the MeV-GeV 7-ray regime. We can expect a similar effect in our theoretical model 
since the escape of the relativistic protons/electrons also occurs at the base of the jet, at the shock 
location. Our model therefore provides a natural connection between the particle acceleration 
mechanism operating in the disk and the formation of the energetic outflows analyzed by Sikora & 
Madejski (2000). 



6. CONCLUSIONS 

In this paper we have demonstrated that particle acceleration at a standing, isothermal shock 
in an ADAF disk can energize the relativistic protons that power the jets emanating from radio- 
loud sources containing black holes. The work presented here represents a new type of synthesis 
that combines the standard model for a transonic ADAF flow with a self-consistent treatment 
of the relativistic particle transport occurring in the disk. The energy lost from the background 
(thermal) gas at the isothermal shock location results in the acceleration of a small fraction of the 
background particles to relativistic energies. One of the major advantages of our coupled, global 
model is that it provides a single, coherent explanation for the disk structure and the formation of 
the outflows based on the well-understood concept of first-order Fermi acceleration at shock waves. 
The theory employs an exact mathematical approach in order to solve simultaneously the coupled 
hydrodynamical and particle transport equations. 

The shock acceleration mechanism analyzed in this paper is effective only in rather tenuous, 
hot disks, and therefore we conclude that our model may help to explain the observational fact 
that the brightest X-ray AGNs do not possess strong outflows, whereas the sources with low X-ray 
luminosities but high levels of radio emission do. We suggest that the gas in the luminous X-ray 
sources is too dense to allow efficient Fermi acceleration of a relativistic particle population, and 
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therefore in these systems, the gas simply heats as it crosses the shock. Conversely, in the tenuous 
ADAF accretion flows studied here, the relativistic particles are able to avoid thermalization due 
to the long collisional mean free path, resulting in the development of a significant nonthermal 
component in the particle distribution which powers the jets and produces the strong radio emission. 

The transport model we have utilized includes the effects of spatial diffusion, first-order Fermi 
acceleration, bulk advection, and vertical escape through the upper and lower surfaces of the disk. 
All of the theoretical parameters can be tied down for a source with a given mass, accretion rate, and 
jet luminosity. The diffusion coefficient model employed in this work is consistent with the behavior 
expected as particles approach the event horizon and also at large radii. Close to the event horizon, 
inward advection at the speed of light dominates over outward diffusion, as expected, while at large 
r, diffusion dominates over advection. The vertical escape of the energetic particles through the 
upper and lower surfaces of the disk at the shock location occurs via spatial diffusion in the tangled 
magnetic field. The approach taken here closely parallels the early studies of cosmic-ray shock 
acceleration. As in those first investigations (e.g., Blandford & Ostriker 1978), we have employed 
the test particle approximation in which the pressure of the accelerated particles is assumed to be 
negligible compared with that of the thermal background gas. This approximation was shown to 
be valid for our specific applications to M87 and Sgr A* in Paper 2. 

We have computed the Green's function describing the particle distribution inside the disk, 
as well as the energy distribution for the relativistic particles that escape at the shock location 
to power the jets. The escaping particle distributions for model 2 (M87) and model 5 (Sgr A*) 
each possess high-energy power-law tails, as expected for a Fermi mechanism. The relatively flat 
slope of the energy distribution also implies the potential for an interesting connection between 
the acceleration process explored here and the formation of the universal cosmic-ray background. 
This is consistent with the results of Szabo & Protheroe (1994), who demonstrated that AGN may 
be an important source of cosmic rays in the region of the "knee." The coupled, self-consistent 
model for the disk structure, the acceleration of the relativistic particles, and the production of the 
outflows presented in this paper represents the first step towards a fully integrated understanding 
of the relationship between accretion dynamics and the formation of the relativistic jets commonly 
observed around radio-loud compact sources. 

The authors would like to thank the referee, Lev Titarchuk, for providing several insightful 
comments that led to significant improvements in the manuscript. 
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APPENDICES 



A. Orthogonality of the Spatial Eigenfunctions 



We can establish the orthogonality of the spatial eigenfunctions Y n (r) by starting with the 
Sturm-Liouville form (see eq. |37| ) 



d 
dr 



dY 



+ \ n co(r)Y n (r) = , 



(Al) 



where ui(r) and S(r) are given by equations (|3Uj) and (|3"5|) . respectively. Let us suppose that X n 
and A m denote two distinct eigenvalues (A„ ^ A m ) with associated spatial eigenfunctions Y n (r) and 
Y m [r). Since Y n and Y m each satisfy equation (|Aip for their respective values of A, we can write 



Yjr) < — 
dr 



S(r) 



w 1 IT, 



S(r) 



dr 
dY r 



dr 



+ X m uj{r) Y m (r) 
+ \ n u(r)Y n (r) 




0. 



Substracting equation ([A3]) from equation (|A2|) yields 



y„(r) £ 

dr 



S(r] 



dY 
(ir 



S(r 



(A n - A m ) w(r) F n (r) F m (r) . 



We can integrate equation (|A4j) by parts from r = r s to r = oo to obtain 



Y n (r)S(r) 



dY 



-Y m (r) S(r] 



dr 



dr 



•s 
oo 



+ 



■ s - ■ s 

■oo 



dY dY 

blr) — — dr 

dr dr 

dY dY 
S( r ) — dr 
dr dr 



(A n - A m ) / w(r) Y n (r) Y m (r) dr . 



(A2) 
(A3) 

(A4) 



(A5) 



Upon cancellation of like terms, this expression reduces to 



\r { \ dY m . s dY n 

Y„(r)— Y m (r 



dr 



dr 



(An -An,) / u{r)Y n {r)Y m {r) dr 



(A6) 



Based upon the asymptotic behaviors of Y n (r) as r — > r s and r — > oo given by equation (|35p . we 
conclude that the left-hand side of equation (jA6j) vanishes, leaving 



u>(r) Y m (r) Y n {r) dr = , m ^ n 



(A7) 



This result establishes that Y m and 1^ are orthogonal eigenfunctions relative to the weight function 
uj{r) defined in equation ([39]) . 



-27- 



REFERENCES 

Abramowicz, M. A. k Chakrabarti, S. K. 1990, ApJ, 350, 281 

Allen S. W., Di Matteo T., Fabian A. C. 2000, MNRAS, 311, 493 

Anyakoha, M. W., Okeke, P. N., k Okoye, S. E. 1987, Ap&SS, 132, 65 

Becker, P. A. 1992, ApJ, 397, 88 

Becker, P. A., & Kafatos, M. 1995, ApJ, 453, 83 

Becker, P. A., k Kazanas, D. 2001, ApJ, 546, 429 

Becker, P. A., & Le, T. 2003, ApJ, 588, 408 

Becker, P. A., k Subramanian, P. 2005 ApJ, 622, 520 

Becker, P. A., Le, T., k Dermer, C. D. 2006, ApJ, 647, 539 

Becker, P. A., Subramanian, P., k Kazanas, D. 2001, ApJ, 552, 209 

Bhattacharyya, S., Bhatt, N., k Misra, R. 2006, MNRAS, 371, 245 

Bhattacharyya, S., Bhatt, N., Misra, R., k Kaul, C. L. 2003, ApJ, 595, 317 

Blandford, R. D., k Begelman, M. C. 1999, MNRAS, 303, LI 

Blandford, R. D., k Ostriker, J. P. 1978, ApJ, 221, L29 

Chakrabarti, S. K. 1989, PASJ, 41, 1145 

Chakrabarti, S. K. 1990, MNRAS, 243, 610 

Chakrabarti, S. K. 1996, ApJ, 464, 664 

Chakrabarti, S. K., Acharyya, K., k Molteni, D. 2004, A&A, 421, 1 
Chakrabarti, S. K., k Das, S. 2004, MNRAS, 349, 649 
Chakrabarti, S. K., & Molteni, D. 1993, ApJ, 417, 671 
Chakrabarti, S. K., k Molteni, D. 1995, MNRAS, 272, 80 
Chen, X., Abramowicz, M. A., k Lasota, J. P. 1997, ApJ, 476, 61 
Das, S., Chattopadhyay, L, k Chakrabarti, S. K. 2001, ApJ, 557, 983 
Dermer, C. D., k Schlickeiser, R. 1993, ApJ, 416, 458 
Di Matteo et al. 2000, MNRAS, 311, 507 



-28- 



Eilek, J. A., & Kafatos, M. 1983, ApJ, 271, 804 

Fermi, E. 1954, ApJ, 119, 1 

Ford, H. C, et al. 1994, ApJ, 435, L27 

Frank, J., King, A. R., & Raine, D. J. 2002, Accretion Power in Astrophysics (Cambridge: Cam- 
bridge University Press) 

Hawley, J. F., Smarr, L. L., & Wilson, J. R. 1984a, ApJ, 277, 296 

Hawley, J. F., Smarr, L. L., k Wilson, J. R. 1984b, ApJS, 55, 211 

Jones, F. C, & Ellison, D. C. 1991, Sp. Sci. Rev., 58, 259 

Kafatos, M., & Yang, R. 1994, MNRAS, 268, 925 

Kazanas, D., Ellison, D. C. 1986, ApJ, 304, 178 

Lanzafame, G., Molteni, D., & Chakrabarti, S. K. 1998, MNRAS, 299, 799 

Le, T., & Becker, P. A. 2004, ApJ, 617, L25 (Paper 1) 

Le, T., & Becker, P. A. 2005, ApJ, 632, 476 (Paper 2) 

Lu, J., Gu, W., & Yuan, F. 1999, ApJ, 523, 340 

Lu, J., & Yuan, F. 1997, PASJ, 49, 525 

Lu, J., & Yuan, F. 1998, MNRAS, 295, 66 

Lynden-Bell, D. 1969, Nature, 223, 690 

Molteni, D., Sponholz, H., & Chakrabarti, S. K. 1996, ApJ, 457, 805 
Morrison, P., Roberts, D., & Sadun, A. 1984, ApJ, 280, 483 
Narayan, R., Kato, S., & Honma, F. 1997, ApJ, 476, 49 

Novikov, I, k Thorne, K. S. 1973, Black Holes, ed. DeWitt, C, & DeWitt, B. S. (New York: New 
York), p. 343 

Owen, F. N., Eilek, J. A., & Kassim, N. E. 2000, ApJ, 543, 611 
Paczyhski, B., & Wiita, P. J. 1980, A&A, 88, 23 
Protheroe, R. J., & Kazanas, D. 1983, ApJ, 265, 620 
Rees, M. J. 1984, Ann. Rev. Astron. Astrophys., 22, 471 



-29- 



Sambruna et al. 2004, ApJ, 608, 698 
Schlickeiser, R. 1989a, ApJ, 336, 243 
Schlickeiser, R. 1989b, ApJ, 336, 264 
Sikora, M., k Madejski, G. 2000, ApJ, 534, 109 
Spruit, H. C. 1987, A&A, 184, 173 

Subramanian, P., Becker, P. A., k Kazanas, D. 1999, ApJ, 523, 203 

Szabo, A. P., k Protheroe, R. J. 1994, Astroparticle Physics, 2, 375 

Titarchuk, L., k Zannias, T. 1998, ApJ, 493, 863 

Urry, C. M., k Padovani, P. 1995, PASP, 107, 715 

Usov, V. V. k Smolsky, M. V. 1998, Ap&SS, 259, 331 

Webb, G. M., Bogdan, T. J. 1987, ApJ, 320, 683 

Yang, R., k Kafatos, M. 1995, A&A, 295, 238 



This preprint was prepared with the AAS IATjrjX macros v5.0. 



